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Abstract. Exact diagonalization is a powerful numerical technique to analyze static 
and dynamical quantities in quantum many-body lattice models. It provides unbiased 
information concerning quantum numbers, energies and wave- functions of the low- 
energy eigenstates for almost any kind of microscopic model. The information about 
energies and quantum numbers is particularly useful to detect possible spontaneous 
symmetry breaking at T = 0. We review some of the advances in the field of frustrated 
quantum magnets which have been possible thanks to detailed symmetry analysis of 
exact diagonalizations spectra. New results concerning the kagome and star lattice 
Heisenberg antiferromagnets are presented. 
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1. Introduction 

There exists a wide range of numerical techniques to deal with quantum many-body 
problems. In the field of lattice spin Hamiltonians, quantum Monte Carlo, density 
matrix renormalization group (DMRG) are among the most powerful, in particular 
because they allow to study accurately very large systems. None of them is however 
efficient to study frustrated models in D > 1. These systems, in D = 2 in particular, 
can present a large variety of phases, with potentially new exotic phases of matter [1]. 

Although limited to small system sizes, exact diagonalizations (ED) give very 
valuable piece of information on these frustrated systems. It can provide (almost) any 
physical quantity: energies, gap, quantum numbers, static and dynamical correlations, 
thermodynamics, ... The method has no bias, and requires no particular knowledge 
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on the low-energy physics. It has been applied from the beginning of numerical 
investigations of quantum spin models.! 

The first question to address in order to characterize a state of matter, concerns the 
(spontaneously) broken symmetry (ies). A natural way to do so, is to evaluate possible 
order parameters from appropriate correlation functions in systems of increasing sizes 
and study if they might remain finite in the thermodynamic limit. Thanks to the rapid 
growth of the computing power, the available system sizes have increased a lot since 
the first ED studies. § Still, these sizes (a few tens of sites at most) are sometimes 
too small to decide if a given order parameter vanishes in the infinite size limit. In 
these cases an approach based on spectrum analysis can be very useful to detect broken 
symmetries. It is simply based on the fact that a spontaneous symmetry breaking (SSB) 
at T = implies some ground-state (quasi-) degeneracy on finite-size spectra. This paper 
illustrates the power of this method through a few examples and references to recent 
studies of 2D quantum antiferromagnets. Some new data on the spin-| Heisenberg 
model on the kagome and "expanded-kagome" lattices (EKL) will also be presented. 

2. Basics of the exact diagonalization method 

This section is a brief presentation of the ED method. It is intended for non specialist 
and is not specific to studies of frustrated magnets. Some readers may therefore wish 
to go directly to Sec. 3. 

2.1. Lanczos method 

The Lanczos method is a numerical iterative algorithm specially efficient to find a few 
extremal eigenvalues and eigenstates (at both ends of the spectrum) for a big sparse 
hermitian matrix. Starting from an initial random vector, one iteratively builds an 
orthonormal basis of the Hilbert space in which the Hamiltonian is tri-diagonal. The 
spectrum of the growing tri-diagonal matrix progressively converges to the spectrum 
of the Hamiltonian. The power of the method comes from the fact that the extremal 
eigenvalues converge first. In particular, the lowest eigenvalue of the tri-diagonal matrix 
converges (exponentially fast) to the true ground-state energy. The number of iterations 
required to obtain the ground-state with a good accuracy is much smaller than the 
Hilbert space dimension. In a typical case, the 2 first eigenvalues of a matrix of size 
10 6 are obtained with machine accuracy in about only 100 iterations (matrix-vector 
multiplications). Once the ground-state has been found, a variant of the Lanczos 
algorithm may be used to compute dynamical correlations (for a recent example, see 
Ref. [5]) or the spin-stiffness [6]. 

| In 1964, Bonner and Fisher [2] studied a by ED 11-site spin-^ chain. 

§ In Ref. [3] the ground-state of a 40-site spin-^ model was computed. A star lattice with 42 sites was 
studied in Ref. [4]. 
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2.2. Other algorithms - Full diagonalization 

Standard algorithms/libraries (like LAPACK) allow to obtain the full spectrum of an 
hermitian matrix. The available Hilbert space dimensions are then however much 
smaller (typically a few 10 4 for a few gigabytes of memory) than with the Lanczos 
method. Full diagonalization allow to compute exactly the finite-temperature properties 
of the system. This may also be a step to solve a generalized diagonalization eigenvalue 
problem as arises when the microscopic spin model has been projected onto a variational 
subspace in a non-orthogonal basis, such as the subspace of first-neighbor valence bond 
coverings (see Refs. [7, 8] and paragraph 3.4 below for an example). 

2.3. Symmetries 

The ED method becomes really powerful only when the spectrum is computed separately 
in each symmetry sector. The symmetry analysis can be described mathematically 
on general grounds. One first determines the symmetry group Q of the Hamiltonian 
H\\ and the irreducible representations (IR) 7o,7i, - "">7n of Q. The projector t\ a 
onto the subspace which transforms according to 7 a (symmetry sector a) reads U a = 
W\ ^g&g Xaig)^ 1 where Xa(g) = T^ilaig)] is the character of 7 Q for the group element 
g (noted g when acting on physical states). By construction, the Il a are orthogonal and 
commute with 7i. Thus one can diagonalize 7i separately in the image of each projector 
(symmetry sectors). In fact, one can further refine the decomposition of the Hilbert 
space. For a given sector a, fl a can be written as a sum: fl Q = J2iZi^a- d a is the 
dimension of the IR 7 Q , the it a are defined by if a = ji ^g^gla'^ig)^ 1 and Ja(g) is 
the i th diagonal element of ^ a {g) (in some arbitrary basis). The ff a also form a family 
of orthogonal projectors commuting with 7i. The spectrum needs not to be computed 
for all subspaces Imf-n^]. One spectrum (e.g. i = 1) will be enough since, in sector a, all 
the 7r^ =1 ' '"' a give the same eigenvalues (hence a degeneracy equal to the IR dimension 
d a )- 

Such a symmetry analysis offers many advantages. The decomposition into stable 
orthogonal subspaces allows to work with smaller vectors and requires less memory and 
CPU time. The density of states being much smaller in each sector than in the full 
spectrum, convergence is reached faster. In the subspace Im[7r^ =1 ], the eigenvalues are 
(generically) non- degenerate. The actual degeneracy of an energy level (hard to find 
with Lanczos if symmetries where not used) is simply given by the IR dimension d a . 
The knowledge of the IR of each energy level is very important to detect SSB. 

In many cases the IR of the full lattice symmetry group can be obtained from one- 
dimensional representations of some subgroups (the IR are induced by some subgroup 
representations). In practice, this amounts first to fix the lattice momentum k and then 
look for IR of its little group (which leaves k invariant). Such a simplified description in 

| In the case of an Heisenberg model on a 2D finite lattice with periodic boundary conditions, this 
group is the direct product of SU(2) (global spin rotations) and the lattice symmetry group. The latter 
usually contains translations as well as point group symmetries (lattice rotations and axis reflections). 
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terms of one-dimensional representations of a subgroup is however not always possible 
(examples: symmetry group of the icosahedron [9].) 

2.3.1. Implementations Usually, a full projection in individual symmetry sectors 
of SU(2) with total spin S is not implemented. This requires large extra CPU 
time/memory and is not necessary. Instead, one merely use a U{1) subgroup, that 
is conservation of S* otai with spin-flip symmetry in the sector S% otgi = which project 
on the subspaces with S even or odd. For SU(2) invariant Hamiltonians, an energy 
level of spin S will then be degenerate in sectors \S Z \ < S which allows to infer 
the value of S. As regards lattice symmetries, a fully automatic treatment can be 
implemented. As an example, the spectra shown in Sec. 3.4 were obtained by a program 
with automatic detection of lattice symmetries and construction of the corresponding 
IR (using GAP[ ]) and application of the appropriate projectors ir. However, to 
our knowledge, no software with such an systematic/automatic treatment of lattice 
symmetries is publicly available at present. 

3. Detecting spontaneous symmetry breaking 

3.1. Quantum numbers 

The existence of a SSB (at T = 0) in the thermodynamic limit, has direct 
consequences on the structure of the low-energy spectrum (degeneracies, quantum 
numbers). Although the SSB only takes place in the thermodynamic limit, this structure 
is often visible on (very) small systems, provided the that the model is not too close to 
a (quantum) phase transition. It is a very useful signature of SSB in ED studies. 

The basic idea is the following. Let be a ground-state of the system and 
Q the symmetry group of the Hamiltonian. If is a broken symmetry state, there 
exists, by definition, at least one group element g G Q under which is not invariant: 
< 1- The linear space V gs generated by all the states {g\ip) , g G Q} has 
thus a dimension d > 1, it defines a (non-trivial) linear representation r of Q. Because 
any g G Q commutes with the Hamiltonian, all the states of V gs are degenerate ground- 
states. The decomposition of Y onto IR 7„ of Q: T = @ a n a ^ ai is obtainable from 
group theory. One finds that the spectrum contains exactly n a ground-state(s) in the 
symmetry sector labeled by the IR a (for n a ^ 0). The multiplicities n a may be obtained 
from the following character representation formula: 



where the states \i) form an orthonormal basis of the ground-state manifold V gs . This 
can also be written as n a = Y2i I |n a |z) 1 1 2 , so that n a > if and only if a broken 
symmetry state has a non-zero projection on sector a. Of course, these multiplicities 
only depend on the symmetry properties of V gs . They can therefore be computed by 
choosing a simpler state, l^o), (without quantum fluctuations for instance) belonging 




(1) 
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to the same "phase" . At this stage the problem of finding the n a has been reduced to a 
purely group-theoretical problem.^ If the system has a finite size, the different sectors 
with n a > will not be exactly degenerate. Still, for a large enough system they should 
become the lowest eigenstates of the spectrum. 

3.2. Discrete broken symmetry, valence-bond crystals 

For a discrete broken symmetry, there are a finite number of degenerate ground-state 
in the thermodynamic limit, separated by a gap (finite in the thermodynamic limit) 
from physical excitations. This (quasi-) degeneracy does not depend on the system 
size (provided boundary conditions do not frustrate the order). In frustrated magnets, 
valence-bond crystal (VBC) are the simplest examples [1]. The procedure described 
in 3.1 gives a straightforward relation between the VBC pattern and the quantum 
numbers of the quasi-degenerate ground-states. For a large enough system, one can 
therefore directly read off the spatial symmetries of the possible VBC from the spectrum, 
without computing any correlation function. To confirm the SSB in the thermodynamic 
limit, one eventually has to check the exponential decay of the energy splitting between 
the quasi-degenerate ground-states. A complementary method based on reduced density 
matrices [ ] can then be used to determine "automatically" the order parameter and 
the VBC pattern. 

Recent ED studies have exhibited various VBC phases. Some 2D examples 
for SU{2) symmetric spin-| models are on the on the honeycomb [12] and square 
lattices [13, 14, 3, 15]. The case of the kagome antiferromagnet will be discussed from a 
VBC point of view in Sec. 3.5 and a new example on the EKL in Sec. 3.4. VBC orders 
may occur in more complex models, but ED remain there limited to quite small sizes. 
For example, ED studies of a 2D SU (4) spin-orbital model on the square lattice indicate 
the possibility of a VBC phase there [16] 

3. 3. Broken continuous symmetry: Neel and nematic phases 

The structure of the spectrum is richer for continuous broken symmetries than for 
discrete ones. First, the existence of gapless Goldstone modes makes the distinction 
between "ground-states" and "excitations" more subtle. And since the ground-state 
subspace (V gs ) has a infinite dimension in the thermodynamic limit (one can perform 
infinitesimal rotations of the initial broken symmetry state), the number of quasi- 
degenerate ground-states has to grow as a function of the system size. 

Frustrated magnets exhibit a remarkably large variety of phases with broken 
continuous symmetries. In a Neel phase, the lattice spontaneously breaks up in 
sublattices in which the all spins point in the same direction (up to quantum zero-point 
fluctuations). The simplest example is the two-sublattice collinear structure realized (at 

% This is particularly simple for one-dimensional representations. In that case one first determines the 
subgroup H of symmetry operations which leaves all the broken symmetry states invariant. Only the 
sectors where Xa(h) = 1 for all h £ H can have n a > 0. 
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T = 0) in the (unfrustrated) antiferromagnetic Heisenberg model on the square lattice. 
But frustration can lead to much more complex structures with non-collinear [17, 3], non- 
coplanar [18] sublattice magnetizations, order-by- disorder effects [19], nematic phases 
with broken SU(2) symmetry but no sublattice magnetization [3, 20], etc. 

To obtain the quantum numbers of the finite-size ground-states in Neel systems, one 
has to decompose a broken-symmetry Neel state \i/jq) (with a well-defined direction for 
the sublattice magnetizations) onto the IR of the SU(2) and lattice symmetry groups. 
As explained above, we can consider a "classical" Neel state with maximum sublattice 
magnetization (each sublattice is fully polarized). For the two-sublattice problem |?/> ) 
can be chosen as an Ising state |^o) — I TIT1 • • •)■ Although trivial, the two-site case 
is instructing: \ip ) = | f|) ~ \S = 0) + \S = 1, S z = 0) is the linear combination 
of a singlet state and a triplet (S = 1) state. Similarly, a general classical Neel state 
for iV (even) sites will have a finite overlap on all total spin sectors from S = to 
S = N/2. Still, the weight of the different spin sectors decreases with increasing S in a 
Neel state (which has no net magnetization in any direction) and one has in particular 
(-00 1 'S'totai = S{S + 1)1^0) ~ N. From this we can guess that, in fact, only spin sectors 
from S = to S ~ \fN are required to construct a state with finite (but not maximum) 
sublattice magnetization and non-zero quantum fluctuations. A simple semi-classical 
argument allows to guess the finite-size scaling of the energies for these states (so- 
called 'Anderson tower" [21], or QDJS in Refs. [17, 22]). As for classical spins, one 
expects a finite uniform susceptibility per site \ m a quantum Neel phase. The total 
magnetization M is therefore given by ~ N\B in presence of an (infinitesimal) applied 
field B. M is also obtained by minimizing E(M) — M ■ B, where E(M) is the energy 
of the lowest state with total spin S — M (in zero field). Combining the two leads to 
E(S) = S 2 /(2Nx)- This spectrum of the Anderson tower corresponds physically to the 
(slow) quantum dynamics of the "rigid body" made by the macroscopic (~ N) sublattice 
magnetizations, it should not be confused with physical spin-wave excitations (higher 
in the spectrum). As expected, the states with S < VN collapse onto the ground-state 
(up to a total energy 0(1)) in the iV — > 00 limit so that a broken symmetry state 
(some combinations of the QDJS) has an energy per site e ~ l/N — > 0. This structure 
(including space group quantum numbers, not discussed here) involves many (~ \/N) 
states and imposes strong constraints on the low energy spectrum. It has been observed 
numerically in finite-size spectra for a large number of quantum antiferromagnets ([2 
and Refs. therein). This tower structure offers a very efficient way to recognize systems 
with continuous broken symmetries, as it already shows up on very small systems (the 
energy gaps in the Anderson tower (~ 1/iV) decay faster than the finite size corrections 
to the sublattice magnetization (~ 1/y/N)). Collinear versus non-collinear Neel states 
can also be readily discriminated from there spectra. A collinear structure has a tower 
with exactly one level per value of S whereas the tower of a non-collinear system (with 
> 3 sublattices) contains 2S + 1 levels in the total spin sector S. We also mention 
that tower structures have also been discussed (and observed numerically) for larger 
symmetry group, like 5'L r (4), in the context of spin-orbital models [23]. 
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A similar analysis can be carried out for quantum p (resp. n) spin-nematics [24], 
which are systems with broken SU(2) symmetry but no sublattice magnetization. A 
simple picture is that the spins spontaneously select a preferred oriented (resp. non- 
oriented) plane, but no particular direction (the order parameter is an anti-symmetric 
(resp. symmetric) rank 2 tensor in the spin variables). Again, such a broken symmetry 
state is a linear combination of many finite-size eigenstates of the systems, with 
specific (and predictable by group theory) total spin and spatial quantum numbers. 
As examples, we mention two recent ED studies which exhibited the tower structures 
associated to a 2D p— [ ] and n— [20] nematics. 

3.4- J e ~ Jt model on the expanded kagome lattice 

The EKL (also dubbed star lattice [ ]) is obtained from the kagome lattice by splitting 
each site into two, and inserting an "expanded" bond connecting the two neighboring 
triangles. Its plaquettes are triangles and dodecagons and all the sites are equivalent. 
Concerning frustration, the classical antiferromagnetic Heisenberg model on the EKL 
has a huge degeneracy equivalent to that on the kagome lattice [4]. Since the EKL has 
two kinds of bonds, (triangle and "expanded" bonds), it is natural to consider a J t — J e 
spin-^ Heisenberg model (both antiferromagnetic here). To our knowledge this model 
has only been studied for J e = J t so far [4]. 

The limit J e ^> J t is simple: the system is made of weakly coupled bonds, the 
ground-state is unique (spin singlet) and all excitations are gapped (A ~ J e ), it is an 
explicit VBC (terminology of Ref. [1]) without any broken symmetry. The model has 
been shown [4] to be in this phase at least up to Jt/J e — 1- 

The limit of weakly coupled triangles (J e <C J t ) is more interesting as the ground- 
state is extensively degenerate at J e = (4 degenerate ground-states per triangle). As 
in Mila's work on the trimerized kagome lattice [25], the degenerate perturbation theory 
for J e /Jt <C 1 can be formulated as a spin-chirality Hamiltonian. In fact, Mila's self- 
consistent mean-field solution of this Hamiltonian almost directly applies on the EKL. 
As a result, one finds an extensive number of (mean-field) ground-states which are in 
one-to-one correspondence with some particular singlet coverings of the EKL (dubbed 
"super coverings" hereafter). These singlet coverings are those which maximize the 
number of occupied J t bonds (as expected since J t 3> J e ). But the extensive degeneracy 
is clearly an artifact of the mean-field approximation. 

To go beyond, we performed some ED of the EKL Heisenberg model restricted the 
subspace of first-neighbor valence bond coverings (a method initiated for the kagome 
antiferromagnet [ ]). As a justification, we note that this variational RVB subspace 
contains both the exact J t = ground-state (a singlet on each J e bond) as well as 
all the mean-field solutions (super coverings) of J e /Jt <C 1. The subspace dimension 
is 2^ Ar+1 ) for N dodecagons and periodic boundary conditions + , which is much smaller 

+ As for the kagome lattice [7, 26], there is a mapping between dimer coverings and Ising pseudo-spin 
configurations. 
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Figure 1. Spectrum of the J e — J< spin-^ Hciscnberg Hamiltonian (on the EKL) 
projected in the first-neighbor RVB subspacc. Left and middle spectra: lattice with 12 
dodecagons (72 sites). Right spectrum: 18 dodecagons (108 sites). The wave-vectors 
symbols are indicated below, in the corresponding Brillouin zones. The VBC pattern 
proposed for J t > 1.3J e is shown on the right of the figure (the singlet "stars", with 
18 sites, are marked with fat bonds). 



than total Hilbert space dimension 2 . Because different singlet coverings are not 
orthogonal, the matrices are not sparse and we had to resort to full diagonalizations. 
Using all lattice symmetries, the complete spectrum has been obtained for N = 12 and 
16 and a part of the spectrum for iV = 18*. The results are summarized Fig. 1. The 
spectrum evolves from a regime with one ground-state and a large gap (Jt/J e ^ 1-1) 
to a regime with many low-energy states (J</J e > 1-3). As expected, the number of 
such low-energy states matches the number of super coverings. For Jt/J e > 1-3 one 
energy level, corresponding to an IR of dimension 2 and momentum in the zone corner, 
is significantly below the other excited states (and goes further down from N = 12 
[middle spectrum] to N = 18 [right spectrum]). This is what one expects for an SSB 
toward a 3-fold degenerate VBC. From the momenta {k = (0,0) and k = (0,±4-7r/3)), 
such VBC should be invariant under any 3-step translations and the simplest guess is a 
crystal with resonating "stars" (analog to that found in the hexagonal lattice quantum 
dimer model [ ]), as shown at the right of Fig. 1. This is also natural because the 
energy gained by singlet resonances is larger of small loops, and the 18-site stars are 
the shortest loops made out (three) super dimers. Although further investigations are 
certainly needed to confirm the existence of this VBC, the spectra (in this variational 
approximation) clearly show that it is the most likely scenario. 

* The two 4-dimensional IR's (subspace dimensions > 29. 10 3 ) corresponding to the 8 wave-vectors 
indicated by 3-leg symbols in the Brillouin zone (right of Fig. 1) could not be computed for N = 18 
(108 sites) because of memory limitations. 
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3.5. Antiferromagnetic spin-^ Heisenberg model on the kagome lattice 

It has been argued that the kagome antiferromagnet could form a VBC at very low 
energies. One scenario [28, 7, 29] is a crystal (VBC-1 in the following), which maximizes 
the number of "perfect" hexagons in a way which preserves the 27r/3 rotation symmetry 
of the lattice (see Fig. 2 of Ref. [ ]). Using the method of Sec. 3.1, we computed the 
quantum numbers associated to this crystal for a 36-site samp left. The lowest exact 
eigenstates (obtained by ED [30]) matching these quantum numbers are marked by * in 
the column VBC-1 of Table I. Although one cannot exclude a deep reorganization of the 
low-energy levels when increasing the system size, these states are not the lowest ones. 
Similarly, a 4- fold degenerate crystal of resonating 12-site "stars" [31] would require the 
two levels indicated in the column "VBC-2" to become the lowest states. The quantum 
numbers of the "columnar" VBC proposed in Ref. [ ] (degeneracy 24) correspond to 
the last column (VBC-3). In all cases, the crystallization would require a complete 
reshuffling of the low-energy spectrum. 
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Table 1. Low energy levels of the Heisenberg model on a 36-site kagome lattice (all are 
total spin singlets) . Momentum k = A refers to the 3 middle points of the (hexagonal) 
zone boundary, B to the two corners of the Brillouin zone. The 6 points C correspond 
to ±B/2. i? 3 (resp. R2) indicates (when applicable) the phase factor acquired by the 
wave-functions under a 27r/3 (resp. it) rotation around the center of an hexagon, a is 
the parity under reflection about the momentum direction. For k — 0, the reflection 
axis tr is chosen parallel to a lattice bond. The eigcnstates which quantum numbers 
match those of VBC ground-states (see text) are marked by * in the last columns. 
There are in fact two kinds of VBC-1, depending on the symmetry of the resonating 
stars (or pin- wheels): even or odd under axis reflections. The 12 levels marked with a 
refer to the "odd-star" VBC-1, and those with b to the "even-star" VBC-1. 
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